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1 Introduction 

Space plasmas are often modeled as a magnctohydrodynamic (MHD) fluid. However, many observed phenomena 
cannot be captured by fluid models, e.g., non-Maxwellian velocity distributions and finite gyro radius effects. 
Therefore kinetic models are used, where also the velocity space is resolved. This leads to a six-dimensional 
problem, making the computational demands of velocity space grids prohibitive. Particle in cell (PIC) methods 
discretize velocity space by representing the charge distribution as discrete particles, and the electromagnetic fields 
are stored on a spatial grid. For the study of global problems in space physics, such as the interaction of a planet 
with the solar wind, it is difficult to resolve the electron spatial and temporal scales. Often a hybrid model is 
then used, where ions are represented as particles, and electrons are modeled as a fluid. Then the ion motions 
govern the spatial and temporal scales of the model. Here we present the mathematical and numerical details of 
a general hybrid model for plasmas. All grid quantities are stored at cell centers on the grid. The most common 
discretization of the fields in PIC solvers is to have the electric and magnetic fields staggered, introduced by Yee 
[17] . This automatically ensures that V • B = 0, down to round-off errors. Here we instead present a cell centered 
discretization of the magnetic field. That the standard cell centered second order stencil for V x E in Faraday's law 
will preserve V • B = was noted by [TJ] . The advantage of a cell centered discretization is ease of implementation, 
and the possibility to use available solvers that only provide for cell centered variables. We also show that the 
proposed method has very good energy conservation for a simple test problem in one-, two-, and three dimensions, 
when compared to a commonly used algorithm. 

2 Definitions 

We have Ni ions at positions r,(£) [m] with velocities Vj(i) [m/s], mass [kg] and charge g, [C], i = 1, . . . , Nj. By 
spatial averaging, we can define the charge density pi(r, t) [Cm -3 ] of the ions, their average velocity u/(r, t) [m/s], 
and the corresponding current density J/(r, f) = piUi [Cm~ 2 s -1 ]. Electrons are modelled as a fluid with charge 
density p e (r,t), average velocity u e (r, i), and current density J e (r,t) = p e u e . The electron number density is 
n e = —p E /e, where e is the elementary charge. If we assume that the electrons are an ideal gas, then p e — n e kT e , 
so the pressure is directly related to temperature (k is Boltzmann's constant). 
The trajectories of the ions are computed from the Lorentz force, 

— =v i5 — = — (E + VixB), t = l,...,Ni 
at at rrii 

where E = E(r, t) is the electric field, and B = B(r,t) is the magnetic field0 



2.1 Hybrid approximations 

A brief overview of hybrid codes can be found in |I6j . A more complete survey can be found in 10 . Most hybrid 
solvers for global simulations have the following assumptions in common. 

1. Quasi-neutrality, pi + p e = 0, so that given the ion charge density, the electron charge density is specified by 
Pe = -pi. 
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1 Usually, charge and current densities are deposited on a grid, using shape functions [3]. 

2 pQ modifies the electric field in the Lorentz force by a term proportional to C and V X B to preserve momentum. 
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2. Ampere's law whithout the transverse displacement current (also called the Darwin approximation, or the 
nonradiative limit) provides the total current, given B, by 



J = [M Q V X B, 

where po = Air ■ 10~ 7 [Hm _1 ] is the magnetic constant (eo/xoc 2 = 1), and from the total current we get 
the electron current, J e = J — J/, and thus the electron velocity, since the quasi- neutrality implies that 

U e = Je/Pe = (J/ - J)/ 'pi- 

3. Massless electrons, m e = 0, lead to the electron momentum equation 

n e m e ^^- = = p e E + J e x B — Vp e + C 
dt 

where the force terms C can be due to collisions, such as electron-ion collisions, electron- neutral [13] collisions, 
or anomalous, i.e. representing electron- wave interactions pQ. In our numerical experiments we have assumed 
that C — I. This provides an equation of state (Ohm's law) for the electric field 

E = — [(J - Jj) X B - Vp e + C] , 
Pi 

with J from Ampere's law. So the electric field is not an unknown. Whenever it is needed, it can be computed. 

4. Faraday's law is used to advance the magnetic field in time, 




= -V x E. 



5. The electron pressure is isotropic (p e is a scalar, not a tensor). 

For the electrons, the remaining degree of freedom is the pressure, p e . Note that p e only affects the ion motions 
through the electric field. The evolution of the magnetic field is not affected since we have V x Vp e = in Faraday's 
law. There are several ways to handle the electron pressure [15l p. 8790], 

1. Assume p e is constant, or zero [5]. 

2. Assume p e is adiabatic (small collision frequency). Then the electron pressure is related to the electron charge 
density by p e cx |p e | 7 , where 7 is the adiabatic index. Commonly used values are 7 = 5/3 [HE], and 7 = 2 

DUG]- 

3. Solve the massless fluid energy equation [TT|l8]. 

^ + u e • Vp e + 7J3e V ■ u e = (7 - 1) ry|J| 2 , 



Here we assume that p e is adiabatic. Then the relative change in electron pressure is related to the relative change 
in electron density by 

Pe_ _ ( "e V 
PeO \neoJ 

where the zero subscript denote reference values. From charge neutrality and p e = n e kT e we have that 

p e = ApJ with A — — p\ 1 T e 

a constant that is evaluated using reference values of pi and T e , e.g., solar wind values. Note that 7 = 1 corresponds 
to assuming that T e is constant, and 7 = gives a constant p e . 
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2.1.1 Hybrid equations 



If we store the magnetic field on a discrete grid Bj, the unknowns are r^, v,, and Bj (supplemented by p e on a grid, 
if we include the electron energy equation) . The time advance of the unknowns can then be written as the ODE 

where Vj x is a discrete rotation operator, and the electric field is 

Ej = — (-J/ x Bj + /i^ 1 (Vj x Bj) x Bj) - Vp e + C. 



3 Discretisation 

An overview of different discretizations of the above equations can be found in 6J Appendix A]. 1 ( Section 3.1] 
provides a consise description of the CAM-CL algorithm introduced by [5]. All our grid variables will be cell 
centered: Bj, Jj, and pj (here Jj and pj are the ionic current and the ionic charge density at cell centers — from 
now we omit the subscript / for simplicity). We follow the Current Advance Method and Cyclic Leapfrog (CAM- 
CL) algorithm [5], but omit the CAM part. The Current Advance Method is used to avoid multiple iterations over 
the particles, but does not conserve energy well, as we will see for a test problem. 

For the particles, we have to solve for and v^; and for the grid cells we have to solve for Bj, Jj and pj. We 
denote time level t = nAt by superscript n. Given B" 1 , r" 1 ' 2 and v™, we do the following steps. 

n+l/2 71-1/2 . \j__n 

n 1 / n+l/2 . n-l/2\ 

r? <~ g ( r i + r * ) ■ 
At r™, deposit particle charges and currents 

Pi->Pj, PiVf^Jj, 

B n+i/2 ^_ B «-i/2 ) p n jn according to CL. (2) 

At r™ +1 ^ 2 , deposit particle charge 

n+l/2 
Pi^Pj 

Estimate electric field at n + 1/2 using the currents at n 

E*^Bf^pfV 2)J n pe) 

vr +1/2 ^v™ + ^^(E* + vrxB; +1 / 2 ). 

At r™ +1 ^ 2 , deposit particle current 

n+l/2 T n+l/2 
^n+1/2 Tjn+1/2 n+l/2 T n+l/2 

v" +1 ^ v? + At— (V l+1/2 + v" +1/2 x B" +1/2N ) . 

Now we have B* +1 , r t t + 1 / 2 anc [ v" +1 . Set n <— n + 1 and start over again. 

For each particle we need a temporary vector. First r™ +1 ^ 2 is temporarily saved during the deposit at r™. Then 
v™ is temporarily saved until the final velocity update. We also need to store the current corresponding to each 
particle, PiV*, in preparation of the deposit operations. 
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The update of the magnetic field in @ using cyclic leapfrog (CL) is done in m sub-time steps of length h = At/m. 
With the notation ~Q P - = Bj ((n + 1/2) At + ph) we have the iteration 

' Bj^Bj-WxEj, 

Bf 1 <- Bf 1 - 2W x E? , p = 1, 2, . . . , m - 1, 

B™ <- B™ _1 — hV x EJ 1 , and 
B „+i/ 2 ^ i ^ B „, + B „\ _ 

Since the magnetic field is leapfrogged in time, we need one temporary grid cell vector. 



3.1 Non-periodic boundary conditions 

To be able to model the interaction of objects with the solar wind, non-periodic boundary conditions in the in- 
direction have been implemented. At x min we have an inflow boundary, and at x max an outflow boundary. The 
other boundaries are still periodic. The computation of V x E in the interior of the simulation domain requires E 
in one extra layer of cells in the ^-directions. Also, computing E in the interior of the simulation domain involves 
V x B, thus also requiring B in one outer layer of cells. At the inflow boundary we specify solar wind values of B 
and E = — uj x B. At the outflow boundary we extrapolate E and B from the interior of the simulation domain 
to one external cell layer (a simple copy of the values from the upstream cells) . 



3.2 Spatial and temporal scales 

If we want solutions of the discrete equations to be accurate approximations of the solutions to the continuous 
equations, a necessary condition is that the discretisation resolves all relevant spatial and temporal scales. The 
smallest spatial scale for the hybrid equations is the ion inertial length (the ion skin depth) 8{ — c/u p i, where c is 
the speed of light and u> P i is the ion plasma frequency, oj pi = riiqf / (eoirii) , n, the ion number density, qi the ion 
charge, rrii the ion mass, and eo ~ 8.854- 10 -12 [Fm _1 ] the vacuum permittivity. The ion inertial length is associated 
with the J x B term in Ohm's law (the Hall term) that describes whistler dynamics. The fastest temporal scale is 
also associated with whistler dynamics. The whistler wave spectrum is cutoff at the electron cyclotron frequency, 
but due to the assumption of massless electrons it is unbounded for the hybrid equations, and the frequency scales 
like uj/Qi = (kc/ujpi) for large k [10] . Here = qiB/uii is the ion gyrofrequency. This gives the CFL constraint 

f Aa 



where n is the spatial dimension. 



4 A quiet plasma test problem 

A uniform, or quiet, plasma is a first test of any simulation code. The solution should only show small statistical 
fluctuations, and energy should be preserved for long simulation times. Matthews [5] describes one- and two- 
dimensional quiet plasma runs, and Brecht [3J present three-dimensional results. 

The number of cells used here is 16, 64 2 , and 32 3 . All boundary conditions are periodic. Ion and electron 
temperatures are given by, = 1, and /3 e = 0. Brecht [3J uses a transport equation for the electron temperature. 
The number of magnetic field sub cycles is 4 in (S], 3 here, and [3j does not use sub cycling. 

Total energy, the sum of the energy stored in the electric and magnetic fields and the kinetic energy of the 
particles, should be conserved. In Table [T] we compare the relative errors in total energy with the published values 
in one-, two-, and three dimensions. 



5 Conclusions 

The hybrid method stores the magnetic field on a grid. Here we have presented a cell centered algorithm as an 
alternative to the staggered grid commonly used. The cell centered method preserves V • B = down to round-off 
errors. In Table [1] it is evident that the proposed method conserves energy well when compared to the commonly 
used CAM-CL method [S|. That the CAM-CL method does not conserve energy well has been noted before [3J[7]. 
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Table 1: Energy errors (total energy) for quiet plasma runs at times T. Numbers in parentheses indicate that the 
parameter was not stated in the reference. 



Reference 


dim. 


particles 
per cell 


Aa; 
St 


At 


T 


error 
Ref. 


error 
Here 




1 


16 


0.5 


0.1 


100 


9% 


0.9% 












300 


47% 


3% 




2 


32 


0.5 


0.1 


100 


2.6% 


0.9% 












300 


14% 


3% 


® 


3 


4 


(1.54) 


0.0056 


112 


<1% 


0.25% 
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